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Abstract 

Using a systematic approach to evaluate Fredholm determinants 
numerically, we provide convincing evidence that the Airyi-process, 
arising as a limit law in stochastic surface growth, is not the limit law 
for the evolution of the largest eigenvalue in GOE matrix diffusion. 

1 Introduction 

One of the unsolved problems in random matrix theory is to understand the 
law for the largest eigenvalue in GOE matrix diffusion. Let M(t) be a matrix- 
valued stationary process on real symmetric matrices of size NxN satisfying: 
(i) the one-time distribution of M(t) is given by the Gaussian Orthogonal 
Ensemble (GOE), (ii) the (independent) entries of M(t) are independent 
stationary Ornstein-Uhlenbeck processes in time. The corresponding process 
for the ordered eigenvalues X^jit), j = 1, . . . , N is Dyson's Brownian motion 
with (3=1. The stationary distribution of the largest eigenvalue, Ajv,jv(*)> 
can be expressed fairly explicit, and its limiting distribution, under proper 
rescaling as iV — > oo, is the GOE Tracy- Widom distribution [16]. However, 
no simple expression for the joint distribution of the largest eigenvalue at 
two different times is known. To be specific, one can ask for the covariance 
of the largest eigenvalue Cov^X^^it), Ajv,jv(0)) and its asymptotic behavior 
as N — > oo. 

For the related case of GUE matrix diffusion, i.e., when M(t) is a hermi- 
tian matrix and the stationary distribution is given by the Gaussian Unitary 
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Ensemble (GUE), the corresponding question is answered. In this case the 
law for the eigenvalues is given by Dyson's Brownian motion with (3 = 2 and 
the limiting process of the properly rescaled largest eigenvalue is the Airy 
process. This process first arose in the context of one-dimensional stochastic 
surface growth with curved macroscopic shape [14] for the so-called polynu- 
clear growth (PNG) model. 

This raised the question, whether the limit process of the largest eigen- 
value in GOE matrix diffusion can also be obtained from a growth process. A 
strong candidate was one-dimensional growth starting from a flat substrate, 
since in this case the limiting one-point distribution is the same as for GOE 
matrix diffusion [13]. This correspondence was partially extended to a mul- 
tilayer version of flat growth with non-intersecting height lines. It was shown 
in [8] that the point process of the multilayer at a fixed position and the 
point process of the GOE ensemble at the edge of the spectrum have the 
same asymptotic law. 

The analogue of the Airy process for flat growth was discovered by 
Sasamoto in a growth model related to PNG [15]. Since its defining ker- 
nel at equal times is, in a certain sense, the square root of the standard Airy 
kernel [10], it was baptized the "Airyi process". Accordingly, for better dis- 
tinction, we call the standard Airy process "Airy 2 process" in the rest of the 
paper. 

In [4] two conjectures have been formulated. The first predicted that the 
Airyx process is also the limit process for the PNG model with flat initial 
conditions, which subsequently has been proven in [6]. The second claimed 
that the Airyi process is also the limit of the largest eigenvalue in GOE 
matrix diffusion ((3=1 Dyson's Brownian motion). 

In this paper we show that the second conjecture does not hold. To 
this end we compare the two-point functions of the Airy! process and of the 
largest eigenvalue in GOE matrix diffusion for different matrix sizes. 

The joint distribution functions for the Airy processes are given in terms 
of Fredholm determinants of integral operators. To evaluate these Fredholm 
determinants we employ a numerical scheme, recently developed by one of 
the authors [2], which in itself is of general interest. For matrix diffusion we 
use straightforward Monte-Carlo simulations on large matrices. 

The comparison shows that the correlation function for GOE matrix dif- 
fusion differs, in the limit of large matrices, from the one for Airyi. In 
contrast, in the case of GUE matrix diffusion, the corresponding numerical 
calculations perfectly illustrate the known convergence to the Airy 2 process. 
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2 Polynuclear growth model 



In this section we present the polynuclear growth model in 1 + 1 dimensions 
and give known results relevant for the discussion. We refer to the original 
papers for more details. 

The model and its multilayer extension 

We briefly define the polynuclear growth (PNG) model on a one-dimensional 
substrate. At time t, the surface is described by an integer-valued height 
function x \— > h(x, t) G Z, x G R, t G R+, with steps of size 1, which is taken 
to be upper semicontinuous, i.e., lim h(x,t) < h(x ,t) for all xo,t. Thus the 



surface consists of up-steps (_T) and down- steps (~L). The dynamics of these 
steps has a deterministic and a stochastic part: 

(i) up- (down-) steps move to the left (right) with unit speed. When a 
down-step and an up-step collide they simply disappear. 

(ii) pairs of up- and down- steps at the same point (spikes) are produced 
by random nucleation events with some given intensity. The up- and 
down-steps of the spikes then spread out with unit speed according to 



The multilayer extension of the PNG model [14] is the following. In- 
stead of a single height function h(x, t) we have a set of height functions 
{hi(x,t),£ < 0}, with the initial condition he(x, 0) = £, for all x G R. The 
dynamics of ho{x, t) is the same as for the original h(x, t). For the remaining 
lines (i) applies as for ho(x,t). Rule (ii) is modified insofar, that for hg(x,t), 
£ < — 1, nucleation events are not produced at random, but whenever there 
is a collision of a pair of steps in level i + 1 at (x, t), a spike is produced in 
level £ at (x, t). By construction the lines do not intersect and one associates 
an (extended) point process i] on R x Z, by 



The PNG droplet 

Consider a flat initial substrate h(x, 0) = 0, x G R. The PNG droplet 
is obtained when the nucleations form a Poisson point process in space- 
time with intensity p(x,t) = 2 for \x\ < t and p(x,t) = otherwise. For 





hg(x, t) = j for some £ < 0, 
otherwise. 



(2.1) 
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large growth time t, the interface has the shape of a droplet, namely the 
deterministic limit, 



/i ma (0 := lim t^hfct, t) = 2 v / l-e%|<i)- (2.2 



t— >oo 



The fluctuations of the height function grow as t 1//3 and the correlation length 
as t 2 / 3 . Therefore, the edge scaling of the (multilayer) height functions around 
the origin, x = 0, is given by 

.droplet, ,x h e (ut 2 / 3 ) - tk^Ut- 1 / 3 ) 

h e («,*):= -fji • (2.3) 

For the PNG droplet, the point process associated to the multilayer is deter- 
minantal. Moreover, rescaled as in (12. 3ft . it converges in the large t limit to 
the Airy field [14], defined by the n-point correlation functions 

p {n \ui, si; . . . ; u n , s n ) = det(K A2 (u u s,; uj, Sj))i<ij<«, (2.4) 

where 



K A2 (u,s; u',s') = < 



poo 

/ dAe ( "'- n)A Ai(s + A)Ai(s' + A), u' < u, 
Jo 

/ clAe (u '-" )A Ai(s + A)Ai(s' + A), v! > u. 

J — oo 



(2.5) 



Denote by A2{u) the highest point of the Airy field at position u. It 
can be seen as a process u t— > A 2 {u) and it is called the Airy 2 process. 
The convergence of the extended point process to the Airy field implies in 
particular that [14] 

lim h* roplct (u,t) = A 2 (u). (2.6) 

The joint distributions of the Airy 2 process are given by Fredholm determi- 
nants: for any given U\ < u 2 < ■ ■ ■ < u m , and s 1 ,...,s m 6R, 

P (^T\{ A ^ U k) < S fc)j = ^i 1 ~ XsK A2 Xs)L\{u 1 ,...,u m }x1R), (2.7) 

where Xs{ u ki x ) = l{ x >s k )- This expression allows to determine some proper- 
ties of the covariance 

g 2 (u) :=Cov(^ a («),^ 2 (0)), (2.8) 

namely 

92 (0) =Vax(A 2 (0)) = 0.81320..., g' 2 (0) = -1, (2.9) 



4 



and the asymptotics for large u [1, 17], 



g 2 (u) = - 2 + - 4 + 0(u~ 6 ), (2.10) 



with the constant c = —3.542 . . . , evaluated numerically from an explicit 
expression in terms of the Hastings-McLeod solution of Painleve II [2]. 

The flat PNG 

Consider a flat initial substrate h(x, 0) = 0, x G R, and run the PNG 
dynamics with constant nucleation intensity, say p(x, t) = 2 for all x G R, 
t > 0. Then the limit shape is flat, /i ma (0 = 2. Thus the edge scaling is 

r M ~*&£=*. (2,1) 

For the flat PNG, the correlation structure of the multilayer version is 
not known, but a few results are available. 

(a) In the large time limit the point process corresponding to (12. lip for 
a fixed value of u converges to a Pfafnan point process [8] whose n-point 
correlation functions are given by 



P 



Mi 



Sl ,...,s n )= 2^ 3 Pf (G GOE (2^ Sl ; 2 2 /%)) 1 < i , i < n . (2.12) 

G is a 2 x 2 matrix kernel (for an explicit expression see, e.g. (2.9) in [8]) 
and Pf is the Pfaffian (Pf(A) = ^det(A) for an antisymmetric matrix A). 
This kernel also occurs for GOE random matrices in the large matrix limit 
at the edge of the spectrum. 

(b) Recently, it has been proven that the Airyi process describes the limit 
of the top line of the multilayer flat PNG [5], namely 

lim (u,t) = 2 1 / 3 A (u/2 2 / 3 ) . (2.13) 

The joint distributions of the Airy! process are given by Fredholm deter- 
minants: for any given u\ < u 2 < ■ ■ ■ < u m , and sj, . . . , s m G R, 

P u\{Ai(u k ) < Sk}^ = det(l - X sK AlX s) L 2 ({iii,...,« m }xR) 1 

(2.14) 

where Xs(uk,x) = t(x > Sk), and the kernel K Al is defined by 
K Al (u,s;u', s') =Ai(s + s'+ (u - u'f) exp ((«' - u)(s + s') + ~{v! - uf) 

' :exp f- ff ~ S)2 ^ \ l(u' > u) (2.15) 



'Airiu' - u) V 4(«'-«) 
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Some properties of the Airy! process like the short and long time behavior 
of the covariance are known. We refer to the review [9] for details. In 
particular, the short-time behavior of the covariance of the Airyi process, 

gi (u) :=Cov(Ai(u),Ai(0)), (2.16) 

satisfies 

^(0) = Var(A(0)) ~ 0.402 . . . , ^(0) = -1. (2.17) 

3 Dyson's Brownian motion 

Dyson's Brownian motion [7] describes the diffusion of N mutually repelling 
particles with positions Xj(t), j = 1, . . . , N, at time t on the real line in a 
harmonic potential, 

dA,-(f) = ( - 7 A,-(f) + § £ A , (t) l A . (t) J dt + db M J = 1. • • • > (3- 1 ) 

the bj(t) being independent standard Brownian motions with Var(6j(£)) = 
t. Let -P(A) denote the probability distribution of particle positions A = 
(Ai, . . . , Ajv)- It satisfies the diffusion equation 

£™ = 1 A; U " I E J^hM - (3.2) 

j=i J \ / J 

The stationary distribution is given by 

P(A) = i|A(A)|^ex P (- 7 f;A J 2 ), (3.3) 

V j= i / 

with A(A) = I j (Xj — Ai), and Z the normalization. 

l<i<j<7V 

In his original work, Dyson linked the special values (5 = 1,2,4 to the 
eigenvalue GOE, GUE and GSE random matrices, respectively. In these 
cases, Xj(t) is the j-th smallest eigenvalue of a random matrix M(t) diffusing 
according to an Ornstein-Uhlenbeck process on real symmetric, hermitian or 
symplectic matrices, respectively. 

We describe the correspondence only in the cases (3=1 (GOE) and 
f3 = 2 (GUE). Let 6g(t), 1 < i,j < N, a = 1,2, be independent standard 
Brownian motions. In the GOE case one sets %(£) = b]j(t) G R, in the 
GUE case one sets b i3 {t) = b}j(t) + ibf^t) E C. Let B^t) = \{kj{t) + 
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bji(t)). Now B(t) is a Brownian motion on the space of real symmetric, 
resp., hermitian matrices, which is invariant with respect to orthogonal, resp. 
unitary rotations. For GOE the independent entries are Bij,j > i, while for 
GUE the independent entries are Ba and Re(Sy), Im(Bij), j > i. These 
real-valued independent entries perform Brownian motions, with variance t 
on the diagonal and variance t/2 for the remaining entries. Now let M(t) be 
the stationary Ornstein-Uhlenbeck process defined by 

dM(t) = - 7 M(t)dt + dB(t). (3.4) 

The stationary distribution is proportional to exp(— 7Tr(M 2 )) in both cases, 
GOE and GUE. By integrating over the angular variables, one gets the 
stochastic evolution of eigenvalues as in (13. 3ft . 

Let us mention here, that the parameter 7 is in fact irrelevant. Multiply- 
ing the eigenvalues by y/j and rescaling time by 7 _1 one can always arrange 
for 7 = 1. We kept this parameter throughout the formulas to facilitate 
comparisons with the literature, where different, sometimes iV-dependent, 
conventions for 7 have been adopted. The most common choice, 7 = 1, leads 
to the standard Hermite kernel with Hermite polynomials orthogonal with 
respect to the weight e~ x . 

GUE diffusion and Airy process 

In the case [3 = 2, the point process associated to the ordered eigenvalues, 
Xj(t) of M(t) is determinantal, defined by the extended Hermite kernel [12]. 
The edge scaling at the upper edge of the spectrum is given by 

\%f(u) = y^iV 1 / 6 (A,( M /( 7 iV 1/3 )) - yfiNh) . (3.5) 

Under this rescaling, the kernel of the corresponding point process converges 
to the Airy kernel K^ 2 as iV — > 00 [1]. This allows to show the convergence 
of the rescaled largest eigenvalue, A^jy (tt), to the Airy 2 process, 

lim \% x %(u) = A 2 (u), (3.6) 

in the sense of convergence of finite-dimensional distributions [11]. The finite- 
iV covariance of the largest eigenvalue is denoted by /^ UE , 

/« UE ( M ) = Cov(A^ E ( M ),A^ E (0)). (3.7) 

Obviously one expects that limAr^oo f^ VE (u) = g2(u). To prove this rigor- 
ously, convergence of moments is needed in (13. 6ft . a result which is currently 
not available. 
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GOE diffusion 



For (3 = 1, the GOE case, explicit expressions for dynamical correlations 
are not known. Nevertheless one expects an analogous behavior as for GUE 
diffusion. In the edge scaling of the ordered eigenvalues Xj(t) of M(t) in the 
GOE case, one has two free parameters (time and space scaling). In order 
to check the hypothesis that the Airy! process describes the evolution of 
the largest eigenvalue for the GOE case, we choose the free parameters by 
matching the covariance and its derivative at zero, see (12.171) . We obtain 

A£2 E H = v^ 176 (A,(2 M /(7iV 1/3 )) - ^Nh) • (3-8) 

As anticipated while speaking of the flat PNG, the point process of the edge- 
rescaled GOE eigenvalues at a fixed time (the one associated to {A^° E (0), 1 < 
j < N}) converges to a Pfaffian point process with n-point correlation given 
by 

pW(si, . . . , s n ) = 2 n Pf (G GOE (2 Sl , 2s j )) 1 < id < n (3.9) 

in the iV — > oo limit. We denote by f GOE the finite- N covariance of the 
largest eigenvalue, 

f% OE (u) = Cov(A^ E ( M ), A^ E (0)). (3.10) 

The scaling in fl£B) is chosen such that f% OE '{0) = -1 and f% OE {0) -> #i(0) 
as N — > oo. As in the GUE case one expects the limit f^ OE (u) = 
limjy^oo f^ OE (u) to exist. In the next section we address the question 
whether f^ OE (u) equals gi(u). 

4 Numerical results 

The Airyi process was regarded as a candidate for the limit of the rescaled 
largest GOE eigenvalue process [4,15], because of the known properties of 
these two processes. Both are stationary processes with the same one-point 
distribution, and the conjectured short time behaviors coincide, too. Fur- 
thermore, as explained above, the underlying multiline ensembles have the 
same limiting single time distribution as point processes on H (see (12.121) and 
(13.91) ). The final ingredient for the guess is that the connections carries over 
the multilines picture in the f3 = 2 case. 

In lack of more analytical input, we looked for an answer to the question, 
whether the Airy! process is the limit of the (3 = 1 Dyson's Brownian motion 
by numerical means. We decided to compare the large N limit of (13. lOf) with 
the covariance of the Airyi process (12.161) . The quantities in question are 
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not straightforwardly accessible. For the Airy! process one needs to evaluate 
Fredholm determinants. For Dyson's Brownian motion we performed Monte- 
Carlo simulations on the eigenvalues of coupled GOE matrices directly. 



Covariances of the Airy! and Airy 2 processes 

The key point of the numerical computation is the evaluation of the Fredholm 
determinants of the two-point joint distributions for the Airy! and Airy 2 pro- 
cesses, eqs. ( 12. 7p . (12.141) . This is explained in details in the recent paper [2]. 
Let us briefly describe the procedure. 

Basic ingredient is a Nystrom-type approximation of integral operators by 
an n-point quadrature formula for integrals over the interval (s, oo) that gives 
exponential convergence rates for holomorphic integrands. Such a formula 
can be based on the holomorphic transformation 

(p s : (0,1) -+ 0,00), £h+ s + 10tan«/2), 

followed by Gauss-Legendre quadrature on the interval (0, 1) with weights 
Wj and points ^ (j = 1, . . . , n). This way we obtain 

f(x) « J^'C(O) f((f) Sk (Q) = ^w kj f{x kj ). 

-k j=l j=l 

The Fredholm determinants (12.71) and (I2.14p are now approximated by the 
mn x mn-dimensional determinant 



det(l-x s ^Xs)i 2 ({ui,..,« ro }xR) ~ det 



(1 -An Ai 2 ••■ A lm \ 
An 1-A 22 ••• A 2m 



\ A m i A m2 • • ■ 1 — A mm J 

(4.1) 

where the submatrices A^ e H nxn (i,j — 1, . . . , m) are defined by 

(Aj)pg = wll 2 K Aii (ui, x ip ; Uj, x jq )w) /2 (p,q=l,..., n). (4.2) 

Theorem 8.1 of [2] shows that the approximation error in (14. ip decays expo- 
nentially with n, that is, like 0(p~ n ) for some constant p > 1. Thus, doubling 
n doubles the number of correct digits; a fact on which simple strategies for 
adaptive error control can be based [3]. Additionally, the level of round-off 
error can be controlled as described in [2]. It turns out that the two-point 
(m = 2) joint distribution can be calculated to an absolute precision of 10~ 14 
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using n quadrature points with n between 20 and 100, depending on the 
specific values of u, si, and 82- The CPU time for a single evaluation of the 
joint distribution is well below a second (using a 2 GHz PC). 

The covariances gi(u) and gi [u] were calculated by first truncating the 
integrals, then integrating by parts (to avoid numerical differentiation), and 
finally using Clenshaw-Curtis quadrature. A truncation at ±10 (except for 
small u < 0.05 in the case of g\, where larger integration intervals are neces- 
sary) and 100 quadrature points in each of the two dimensions are sufficient 
to secure an absolute precision of 10~ 10 . 

The code for calculating the two-point joint distributions and the covari- 
ances gi and gi can be obtained from the first author upon request. 

Monte Carlo Simulation of Random Matrices 

To get an estimate for the covariance of the largest eigenvalue of GUE and 
GOE matrices, we performed straightforward Monte-Carlo simulations. A 
collection of matrices C*., k — 0, . . . , K, independently distributed according 
to the stationary distribution of (13.41) can be easily produced with standard 
pseudo random generators, since each Ck consists of independently normally 
distributed entries. Fixing a time step At, it is easy to see that for the station- 
ary process M(t) governed by (13. 4j) the joint distribution of (M(kAt)) 0<k<K 
is the same as for the matrices M^, defined by 

M = C , M k = e-^Mfc.! + \A - e-^ A 'C fc , 1 < k < K. (4.3) 

Now, one numerically determines the largest eigenvalues of the and 
rescales according to eqs. (I3.5p . resp., (13. 8p . This yields realizations of 
^n°n( u ) an d ^ < n J n( u ) a ^ equidistant times u. The empirical auto-covariance 
of these time series gives an estimate for the covariance functions f^ OE (u) 
and /^ ue (m) at discrete values of u. The data presented here are for N = 64, 
128, and 256 with 7 = 1/2 and At = ^N' 1 ^. We chose K = 10 6 , and 
produced up to 100 independent realizations of the time series in each case. 
This allows to obtain an error estimate for each data point. 

Discussion 

In Figure [1] we compare the covariance (12.161) of the Airyi process and the 
one of the largest eigenvalue for GOE matrix diffusion (13. lOj) for iV = 64 and 
iV = 256. One clearly sees that they do not agree. 

Increasing the matrix dimension does not change sensibly the result; 
namely, the results for N = 128 agree to plotting accuracy with the one 
for iV = 256 in Figure [1] and therefore are not shown. In comparison, in 
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GOE versus Airvj 



0.45 




u 



Figure 1: The covariance g\ of the Airyi process (line) versus the one of the 
largest eigenvalue for GOE matrix diffusion, f^ OE , for N = 64, 256. 



Figure [2] we plot the covariance ( 12. 8p of the Airy2 process and the one of 
the largest eigenvalue for GUE matrix diffusion ( 13. 7p for the same matrix di- 
mensions. Here the agreement is already quite good even for relatively small 
matrix sizes. In both plots the errorbars are of order 10~ 3 , smaller then the 
symbols used and therefore omitted. 

Concerning the decay of g\, it appears to be superexponentially in sharp 
contrast to the algebraic decay ( 12.101) of g%. After reinspecting the 2x2 
block Fredholm determinant this behavior becomes clear, since one of the 
off-diagonal blocks is superexponentially small in u for large values of it, 
while the others stay of order one, a fact already noticed by Widom [18]. 

Finally, in Figure El we provide a comparison of the decay of correlation 
for GOE and GUE matrix diffusion. In a log-log plot we draw and 
f^ VE for iV = 128 and N = 256 with errorbars. For GUE one observes 
the deviation from the asymptotic behavior u~ 2 for large u due to finite size 
effects. Remarkably f^ OE (u) looks very similar to ^f^ VE (2u), indicating 
that the large u behavior of f^ OE {u) might also be algebraically decaying, 
in sharp contrast to the superexponential decay of gi(u). Given the small 
matrix dimensions we used, we can not, however, conclude whether the decay 
for f^ OE {u) is of order u~ 2 as for GUE or not. 
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GUE versus Airy 2 




1 2 3 4 5 



u 

Figure 2: The covariance g-i of the Airy 2 process (line) versus the one of the 
largest eigenvalue for GUE matrix diffusion, /^ UE , for N = 64, 256. 



GOE versus GUE 
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Figure 3: Log-log plot of the rescaled correlation functions for GOE and 
GUE. 
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